library(sfsmisc)
setwd("C:\\Users\\Ofer\\Desktop\\My Dropbox\\PhD\\LATEX\\4\\figures")
#par(mfrow=c(2,1))
alpha<-.5;min<-1;max<-100;n<-10000;col="black";add=F
#alpha<-.5;min<-1;max<-100;n<-10000;col="blue";add=T
#alpha<-1;min<-1;max<-100;n<-10000;col="red";add=T
#alpha<-5;min<-1;max<-100;n<-10000;col="green";add=T

####################
# plot function 
f.dist<-function(x, min=1, max=10, alpha=1.4){alpha*x^(-alpha-1)/(min^-alpha - max^-alpha)}
png("ImgPowerDist.png",  bg = "transparent")
# exact function of the distribution
curve(f.dist(x, min=min, max=max, alpha=alpha), 
	from=min, to=max, col=col, lwd=2,
    	add=add,
	xlim=c(0,100),
	ylab='',  lty = 2, 
	xlab='', pch=20, cex=.1,  xaxt = "n",  yaxt = "n")
mtext("frequency", side=2, line=2.5, cex=1.5)
mtext("x-value", side=1, line=2.5, cex=1.5)
eaxis(1, cex.axis=.8)
eaxis(2, cex.axis=.8)

# integrate(f, lower=min, upper=max, min=min, max=max, alpha=alpha)
# plot distribution
generatePower<- function(nmax, alpha, min=1, max=Inf) {
	u<-runif(nmax) # alpha<-1; min=1; max=10; u<-.001
	norm<-(min^-alpha - max^-alpha)
	(min^-alpha - u*norm)^-(1/alpha)
}
d<-generatePower(n, alpha, min=min, max=max)
hist<-hist(d, plot = FALSE, breaks=100)
points( hist$mids, hist$counts/n, 
	col="black", pch=3, cex=1, ylab='', xlab='x value', 
	# log="xy",
	xaxt = "n",  yaxt = "n")

# plot cumulative distribution
f.cum<-function(x, min=1, max=10, alpha=1.4){alpha*(x^-alpha-max^-alpha)/(min^-alpha - max^-alpha)}
curve(2*f.cum(x, min=min, max=max-1e-9, alpha=alpha), 
	from=min, to=max, col="red", lwd=2, add=T)

normalized<-rev(cumsum(rev(hist$counts/n))) # /  sum(rev(cumsum(rev(hist$counts/10000))))
points( hist$mids, normalized  , 
 	col="blue", pch=4, cex=1, ylab='', xlab='x value')

#legend
leg.txt <- c("Exact function (DF)",
		 "Simulation     (DF)", 
		 "Exact function (reverse CDF)", 
             "Simulation     (reverse CDF)")
legend(x=15, y=.5, leg.txt,  lty = c(2, -1, 1, -1), pch = c(-1, 3, -1, 4),
      col = c("black", "black", "red", "blue"), 
	cex = 1, lwd=c(2,1,2,1))

dev.off()

#############################################################


#############################################################
# plot log cumulative distribution
library(sfsmisc)
setwd("C:\\Users\\Ofer\\Desktop\\My Dropbox\\PhD\\LATEX\\4\\figures")
alpha<-.5;min<-1;max<-100;n<-10000;col="black";add=F
f.dist<-function(x, min=1, max=10, alpha=1.4){alpha*x^(-alpha-1)/(min^-alpha - max^-alpha)}
f.cum<-function(x, min=1, max=10, alpha=1.4){(x^-alpha-max^-alpha)/(min^-alpha - max^-alpha)}
png("ImgPowerDistLogLog.png",  bg = "transparent")

curve(f.cum(x, min=min, max=max-1e-9, alpha=alpha), 
	from=min, to=max, col="red",
    	# add=T,
	ylim=c(1e-4, 1), lwd=2,
	ylab='', xlab='', 
	pch=20, cex=.01,  xaxt = "n",  yaxt = "n" , log="xy")
mtext("frequency", side=2, line=2.5, cex=1.5)
mtext("x-value", side=1, line=2.5, cex=1.5)
eaxis(1, cex.axis=.8)
eaxis(2, cex.axis=.8)


cumNormalized<-rev(cumsum(rev(hist$counts/n))) 
points( hist$mids, cumNormalized, 
	col="blue", pch=4, cex=1, ylab='', xlab='x value', xaxt = "n",  yaxt = "n")
#eaxis(1)  axTicks(1)[c(T, F, F, T, F, F, T)])
#eaxis(2, axTicks(2)[c(T, T, T, T)])

# plot log distribution
points(hist$mids, hist$counts/n,  
	# xaxt = "n",  yaxt = "n", log="xy", 
	col="black", pch=3, cex=1)

curve(f.dist(x, min=min, max=max-1e-10, alpha=alpha), col="black", lty=2,
	pch=20, cex=.01, add=T, lwd=2)
#legend
leg.txt <- c("Exact function (DF)",
		 "Simulation     (DF)", 
		"Exact function (reverse CDF)", 
             "Simulation     (reverse CDF)" )
legend(x=1, y=1e-3, leg.txt, lty = c(2, -1, 1, -1), pch = c(-1, 3, -1, 4),
      col = c("black", "black", "red", "blue"), 
	cex = 1, lwd=c(2,1,2,1))

dev.off()

#==============================================================0

f.cum<-function(x, min=1, max=10, alpha=1.4)
{(x^-alpha-max^-alpha)/(min^-alpha - max^-alpha)}


f.log.approx<-function(x, min=1, max=10, alpha=1.4){alpha*log(min)-alpha*log(x)}

{
	r<-(min/max)^alpha; 
	y<-(1+(x-min)/min)^alpha
	# y<-1+alpha*(x-min)/min
	
	list(y0, y, ret)
}


{( (max^alpha-x^alpha)*(min*max)^alpha) / ((max^alpha-min^alpha)*(x*max)^alpha) }
{( (1-(x/max)^alpha)*(min/max)^alpha) / ((1-(min/max)^alpha)*(x/max)^alpha) }

x<-min+1
f.cum(x,min, max, alpha)
f.approx(x,min, max, alpha)

f.log.approx(x,min, max, alpha)
log(f.cum(x,min, max, alpha))

